Libraries

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──
## ✔ ggplot2 3.4.1     ✔ purrr   1.0.1
## ✔ tibble  3.1.8     ✔ stringr 1.5.0
## ✔ tidyr   1.3.0     ✔ forcats 1.0.0
## ✔ readr   2.1.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(sp)
library(ggplot2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(elevatr)
library(ggspatial)
library(ggmap) #citation("ggmap")
## ℹ Google's Terms of Service: <]8;;https://mapsplatform.google.comhttps://mapsplatform.google.com]8;;>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(leaflet)
library(htmlwidgets)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(maps)
## 
## Attaching package: 'maps'
## 
## The following object is masked from 'package:purrr':
## 
##     map
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(ggsn)
## Loading required package: grid
## 
## Attaching package: 'ggsn'
## 
## The following object is masked from 'package:raster':
## 
##     scalebar
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
library(performance)
library(ggResidpanel)
library(ape)
## Registered S3 method overwritten by 'ape':
##   method   from 
##   plot.mst spdep
## 
## Attaching package: 'ape'
## 
## The following objects are masked from 'package:raster':
## 
##     rotate, zoom
## 
## The following object is masked from 'package:dplyr':
## 
##     where
library(reshape2)
## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths

Data

# Colors to Use
ryg.palette<-c("#DB4325", "#EDA247", "#E6E1BC","#57C4AD", "#006164")
ryg.palette.long<-c("#DB4325", "#EDA247", "#E6E1BC","#57C4AD", "#006164", "grey")
ryg.palette.short<-c("#EDA247", "#E6E1BC","#57C4AD", "#006164")

# API Key for Geocoding:
api_key = "AIzaSyAeoVaM_45SQ4G-lHeci1RhRNbhpxO3BDc"
register_google(key = api_key) 

waveone<-read.csv("/Volumes/research/Marshall Fire Health/marshall_w1_deID.csv")
wavetwo<-read.csv("/Volumes/research/Marshall Fire Health/marshall_w2_deID.csv")
wave.one<-waveone %>%
  mutate(mailingcity=recode_factor(mailingcity, "Superior"="SUPERIOR", "UNINCORPORATED"="BOULDER")) %>%
  mutate(across(air_quality_1:air_quality_4, as.factor)) %>%
  mutate(across(remediation_1:remediation_4, as.factor)) %>%
  mutate(across(air_cleaning_1:air_cleaning_4, as.factor)) %>%
  mutate(group=as.factor(group)) %>% 
  mutate(ownership_status=recode_factor(ownership_status,
                                        "1"="Owned",
                                        "2"="Rented",
                                        "3"="Owned but Not Living",
                                        "4"="Purchased after Dec. 30, 2021",
                                        "5"="Other")) %>%
  mutate(air_quality_1=recode_factor(air_quality_1,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>%
  mutate(air_quality_2=recode_factor(air_quality_2,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>%
  mutate(air_quality_3=recode_factor(air_quality_3,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>%
  mutate(air_quality_4=recode_factor(air_quality_4,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>% 
  mutate(home_smell_1week=recode_factor(home_smell_1week,
                                        "0"="No",
                                        "1"="Yes",
                                        "2"="Not sure/don't remember")) %>%
  mutate(home_smell_type=recode_factor(home_smell_type,
                                       "1"="Like a campfire",
                                       "2"="Like chemicals or a chemical fire",
                                       "3"="Other")) %>%
  mutate(group=recode_factor(group,
                             "A"="within fire perimeter",
                             "B"="within 1/2 mile",
                             "C"="1/2 to 1 mile",
                             "D"="1 to 2 miles")) %>%
  mutate(remediation_1=recode_factor(remediation_1, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither"))%>%
  mutate(remediation_2=recode_factor(remediation_2, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither"))%>%
  mutate(remediation_3=recode_factor(remediation_3, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither"))%>%
  mutate(remediation_4=recode_factor(remediation_4, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither"))%>%
  mutate(air_cleaning=recode_factor(air_cleaning_1, 
                                    "1"="Changing air filters in HVAC")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_2, 
                                    "1"="Using air cleaners"))%>%
  mutate(air_cleaning=recode_factor(air_cleaning_3, 
                                    "1"="Hiring someone to clean indoor air")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_4, 
                                    "1"="Other")) %>%
  mutate(air_cleaning_3=recode_factor(air_cleaning_3,
                                      "1"="Hiring someone to clean indoor air"))%>%
  mutate(air_cleaning_4=recode_factor(air_cleaning_4,
                                      "1"="Other")) %>%
  mutate(education=recode_factor(education,
                                 "0"="Did not specify", 
                                 "1"="High School Graduate",
                                 "2"="GED or Equivalent",
                                 "3"="Some College",
                                 "4"="Associate's Degree",
                                 "5"="Bachelor's Degree",
                                 "6"="Master's Degree",
                                 "7"="Doctoral Degree")) %>%
  mutate(income=recode_factor(income,
                              "1"="Less than $20,000", 
                              "2"="$20,000 to $34,999",
                              "3"="$35,000 to $49,999",
                              "4"="$50,000 to $79,999",
                              "5"="$80,000 to $99,999",
                              "6"="$100,000 to $149,999",
                              "7"="$150,000 to $199,999",
                              "8"="$200,000 or more",
                              "9"="Prefer not to answer")) %>%
  filter(finished==1) %>% #only analyzing finished data, not in wave.two
  filter(mailingcity=="BOULDER"|mailingcity=="SUPERIOR"|mailingcity=="LOUISVILLE") #filtering only respondents in Boulder Unincorporated, Superior or Louisville

wave.two<-wavetwo %>%
  mutate(air_quality_2=as.factor(air_quality_2)) %>%
  mutate(air_quality_4=as.factor(air_quality_4)) %>%
  mutate(across(remediation_1:remediation_4, as.factor)) %>%
  mutate(across(air_cleaning_1:air_cleaning_4, as.factor)) %>%
  mutate(group=as.factor(group)) %>% 
  mutate(renterowner=recode_factor(renterowner,
                                        "1"="Owned",
                                        "2"="Rented")) %>%
  mutate(air_quality_2=recode_factor(air_quality_2,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>%
  mutate(air_quality_4=recode_factor(air_quality_4,
                                     "1"="Strongly disagree",
                                     "2"="Somewhat disagree",
                                     "3"="Neither agree nor disagree",
                                     "4"="Somewhat agree",
                                     "5"="Strongly agree",
                                     "6"="Don't know")) %>%
  mutate(group=recode_factor(group,
                             "A"="within fire perimeter",
                             "B"="within 1/2 mile",
                             "C"="1/2 to 1 mile",
                             "D"="1 to 2 miles")) %>%
  mutate(remediation_1=recode_factor(remediation_1, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither")) %>%
  mutate(remediation_2=recode_factor(remediation_2, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither")) %>% 
  mutate(remediation_3=recode_factor(remediation_3, 
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither")) %>%
  mutate(remediation_4=recode_factor(remediation_4,
                                     "1"="Did this already",
                                     "2"="Plan to do this",
                                     "3"="Neither")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_1, 
                                    "1"="Changing air filters in HVAC")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_2, 
                                    "1"="Using air cleaners")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_3, 
                                    "1"="Hiring someone to clean indoor air")) %>%
  mutate(air_cleaning=recode_factor(air_cleaning_4,
                                    "1"="Other")) %>%
  mutate(income=recode_factor(income,
                              "1"="Less than $20,000", 
                              "2"="$20,000 to $34,999",
                              "3"="$35,000 to $49,999",
                              "4"="$50,000 to $79,999",
                              "5"="$80,000 to $99,999",
                              "6"="$100,000 to $149,999",
                              "7"="$150,000 to $199,999",
                              "8"="$200,000 or more",
                              "9"="Prefer not to answer")) %>%
  mutate(education=recode_factor(education,
                                 "0"="Did not specify", 
                                 "1"="High School Graduate",
                                 "2"="GED or Equivalent",
                                 "3"="Some College",
                                 "4"="Associate's Degree",
                                 "5"="Bachelor's Degree",
                                 "6"="Master's Degree",
                                 "7"="Doctoral Degree"))

Descriptive Survey Statistics (tables)

Number of respondents for Wave 1 and Wave 2

nrow(waveone)
## [1] 3462
nrow(wavetwo)
## [1] 576

Number of respondants who completed the survey for Wave One and Wave Two

Respondants = c(nrow(wave.one),nrow(wave.two))
Wave=c("One","Two")
df1 = cbind(Wave,Respondants) %>% 
  as.data.frame()
knitr::kable(df1)
Wave Respondants
One 802
Two 576

% male/female

Wave One

tblFun(wave.one$Gender3)
Count Percentage
Female 473 58.98
Male 298 37.16
Other or Declined 31 3.87

Wave Two

tblFun(wave.two$Gender3)
Count Percentage
Female 351 60.94
Male 201 34.90
Other or Declined 24 4.17

Educational attainment for each wave

Wave One

tblFun(wave.one$education)
Count Percentage
High School Graduate 13 1.65
GED or Equivalent 1 0.13
Some College 43 5.46
Associate’s Degree 21 2.66
Bachelor’s Degree 310 39.34
Master’s Degree 283 35.91
Doctoral Degree 117 14.85

Wave Two

tblFun(wave.two$education)
Count Percentage
High School Graduate 9 1.59
GED or Equivalent 1 0.18
Some College 28 4.96
Associate’s Degree 13 2.30
Bachelor’s Degree 213 37.70
Master’s Degree 216 38.23
Doctoral Degree 85 15.04

Income for each wave

Wave One

tblFun(wave.one$income)
Count Percentage
Less than $20,000 7 0.89
$20,000 to $34,999 11 1.40
$35,000 to $49,999 27 3.45
$50,000 to $79,999 50 6.39
$80,000 to $99,999 62 7.92
$100,000 to $149,999 148 18.90
$150,000 to $199,999 125 15.96
$200,000 or more 210 26.82
Prefer not to answer 143 18.26

Wave Two

tblFun(wave.two$income)
Count Percentage
Less than $20,000 5 0.89
$20,000 to $34,999 9 1.60
$35,000 to $49,999 19 3.37
$50,000 to $79,999 32 5.67
$80,000 to $99,999 45 7.98
$100,000 to $149,999 110 19.50
$150,000 to $199,999 98 17.38
$200,000 or more 150 26.60
Prefer not to answer 96 17.02

Race/ethnicity

Wave One

tblFun(wave.one$RaceEthn2)
Count Percentage
35 4.36
POC 76 9.48
White 691 86.16

Wave Two

tblFun(wave.two$RaceEthn2)
Count Percentage
23 3.99
POC 46 7.99
White 507 88.02

Owner/Renter

Wave One

tblFun(wave.one$ownership_status)
Count Percentage
Owned 744 93.12
Rented 31 3.88
Owned but Not Living 10 1.25
Purchased after Dec. 30, 2021 4 0.50
Other 10 1.25

Wave Two

tblFun(wave.two$renterowner)
Count Percentage
Owned 554 96.18
Rented 22 3.82

Impact category

Wave One

tblFun(wave.one$impact_cat)
Count Percentage
7 0.87
Complete loss 204 25.44
Damaged, living there 355 44.26
Damaged, not living there 41 5.11
No damage, living there 191 23.82
No damage, not living there 4 0.50

Wave Two

tblFun(wave.two$impact_cat)
Count Percentage
Complete loss 163 28.30
Damaged, living there 249 43.23
Damaged, not living there 33 5.73
No damage, living there 128 22.22
No damage, not living there 3 0.52

Distance category

Wave One

tblFun(wave.one$group)
Count Percentage
within fire perimeter 465 57.98
within 1/2 mile 197 24.56
1/2 to 1 mile 74 9.23
1 to 2 miles 66 8.23
0 0.00

Wave Two

tblFun(wave.two$group)
Count Percentage
within fire perimeter 343 59.55
within 1/2 mile 139 24.13
1/2 to 1 mile 56 9.72
1 to 2 miles 38 6.60

Air Quality Perceptions

Wave One:

Wave One only: compare before and after the fire – perception of air quality in their neighborhood

#trying a vertical graph
w1.before.neigh.plot=ggplot(filter(wave.one,!is.na(air_quality_1)), aes(fill=air_quality_1,x="Before the Fire")) +
  geom_bar(position="stack") +
  labs(title = "I am confident that the air in my neighborhood is safe to breathe") +
  xlab(NULL) +
  ylab(NULL) +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

w1.after.neigh.plot=ggplot(filter(wave.one,!is.na(air_quality_3)), aes(fill=air_quality_3,x="After the Fire")) +
  geom_bar(position="stack") +
  xlab(NULL) +
  ylab("Count") +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

grid.arrange(w1.before.neigh.plot,w1.after.neigh.plot)

Wave One only: compare before and after the fire – perception of air quality in their neighborhood by distance

ggplot(filter(wave.one,!is.na(air_quality_1)), aes(fill=air_quality_1, x=group)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.one,!is.na(air_quality_3)), aes(fill=air_quality_3, x=group)) +
    geom_bar(position="stack") +
  labs(title = "Currently, I am confident that the air in my neighborhood is \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One only: compare before and after the fire – perception of air quality in their neighborhood by impact category

ggplot(filter(wave.one,!is.na(air_quality_1)), aes(fill=air_quality_1, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.one,!is.na(air_quality_3)), aes(fill=air_quality_3, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Currently, I am confident that the air in my neighborhood is \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One only: compare before and after the fire – perception of air quality in their home

#trying a horizontal graph
w1.before.home.plot=ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2,y="Before the Fire")) +
  geom_bar(position="stack") +
  labs(title = "I am confident that the air inside my home is safe to breathe") +
  ylab(NULL) +
  xlab(NULL) +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

w1.after.home.plot=ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4,y="After the Fire")) +
  geom_bar(position="stack") +
  ylab(NULL) +
  xlab("Count") +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

grid.arrange(w1.before.home.plot,w1.after.home.plot)

Wave One only: compare before and after the fire – perception of air quality in their home by distance

ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2, x=group)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4, x=group)) +
    geom_bar(position="stack") +
  labs(title = "Currently, I am confident that the air in my neighborhood is \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One only: compare before and after the fire – perception of air quality in their home by impact category

ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Currently, I am confident that the air in my neighborhood is \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One and Two:

Wave One to Wave Two comparison – perception of air quality in neighborhood (AQ2)

w1.after.neigh.plot=ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2,y="Wave One")) +
  geom_bar(position="stack") +
  labs(title = "I am confident that the air in my neighborhood is safe to breathe") +
  ylab(NULL) +
  xlab(NULL) +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

w2.after.neigh.plot=ggplot(filter(wave.two,!is.na(air_quality_2)), aes(fill=air_quality_2,y="Wave Two")) +
  geom_bar(position="stack") +
  ylab(NULL) +
  xlab("Count") +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

grid.arrange(w1.after.neigh.plot,w2.after.neigh.plot)

Wave One to Wave Two comparison – perception of air quality in neighborhood (AQ2) by distance

w1.after.neigh.plot=ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2,y="Wave One")) +
  geom_bar(position="stack") +
  labs(title = "I am confident that the air in my neighborhood is safe to breathe") +
  ylab(NULL) +
  xlab(NULL) +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

w2.after.neigh.plot=ggplot(filter(wave.two,!is.na(air_quality_2)), aes(fill=air_quality_2,y="Wave Two")) +
  geom_bar(position="stack") +
  ylab(NULL) +
  xlab("Count") +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

grid.arrange(w1.after.neigh.plot,w2.after.neigh.plot)

Wave One to Wave Two comparison – perception of air quality in neighborhood (AQ2) by impact category

ggplot(filter(wave.one,!is.na(air_quality_2)), aes(fill=air_quality_2, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.two,!is.na(air_quality_2)), aes(fill=air_quality_2, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "Before the Marshall Fire, I was confident that the air in my neighborhood was \n safe to breathe",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One to Wave Two comparison – perception of air quality in home (AQ4)

w1.after.home.plot=ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4,y="Wave One")) +
  geom_bar(position="stack") +
  labs(title = "I am confident that the air in my neighborhood is safe to breathe") +
  ylab(NULL) +
  xlab(NULL) +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

w2.after.home.plot=ggplot(filter(wave.two,!is.na(air_quality_4)), aes(fill=air_quality_4,y="Wave Two")) +
  geom_bar(position="stack") +
  ylab(NULL) +
  xlab("Count") +
  scale_fill_manual(values = ryg.palette) +
  theme(plot.title = element_text(hjust = 0.5),   
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill=guide_legend(title=NULL))

grid.arrange(w1.after.neigh.plot,w2.after.neigh.plot)

Wave One to Wave Two comparison – perception of air quality in home (AQ4) by distance

ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4, x=group)) +
    geom_bar(position="stack") +
  labs(title = "I am confident that the air inside my home was \n safe to breathe",
       subtitle = "Wave One",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.two,!is.na(air_quality_4)), aes(fill=air_quality_4, x=group)) +
    geom_bar(position="stack") +
  labs(title = "I am confident that the air inside my home was \n safe to breathe",
       subtitle = "Wave Two",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Wave One to Wave Two comparison – perception of air quality in home (AQ4) by impact category

ggplot(filter(wave.one,!is.na(air_quality_4)), aes(fill=air_quality_4, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "I am confident that the air inside my home was \n safe to breathe",
       subtitle = "Wave One",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

ggplot(filter(wave.two,!is.na(air_quality_4)), aes(fill=air_quality_4, x=impact_cat)) +
    geom_bar(position="stack") +
  labs(title = "I am confident that the air inside my home was \n safe to breathe",
       subtitle = "Wave Two",
       x = "Distance from fire perimeter", fill ="") +
   scale_fill_manual(values = ryg.palette)

Maps:

Maps Data:

Map of perception of air quality in one’s neighborhood - color from stacked plot colors

Before the Fire

geo_marshall_1=filter(geo_marshall,!is.na(air_quality_1)) %>% 
  st_as_sf(coords=c("lon","lat"), crs = 4326)
marshall.sp.aq1 = as(geo_marshall_1,"Spatial")
factpal <- colorFactor(c("#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641"), c("Strongly disagree","Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"))

leaflet(marshall.sp.aq1) %>%
  addTiles() %>%
  addCircleMarkers(group=marshall.sp.aq1$air_quality_1, color = ~factpal(air_quality_1),
                   radius = 1,
                   opacity = 1,
                   lng = marshall.sp.aq1@coords[,1],
                   lat = marshall.sp.aq1@coords[,2])
ggplot(filter(geo_marshall_1,!is.na(air_quality_1))) +
  annotation_map_tile() +
  annotation_scale() +
  geom_sf(aes(color=as.factor(air_quality_1))) +
  scale_color_manual(values = c("Strongly disagree"="#d7191c",
                                "Somewhat disagree"="#fdae61",
                                "Neither agree nor disagree"="#ffffbf",
                                "Somewhat agree"= "#a6d96a",
                                "Strongly agree" = "#1a9641")) +
  labs(title="Marshall Fire Survey Results (NA's filtered out)", 
       subtitle=paste0("Wave One: Total Respondants Displayed: ",nrow(filter(geo_marshall_1,!is.na(air_quality_1)))),
       caption="Statement: Before the Marshall Fire, I was confident that the \n air in my neighborhood was safe to breath") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5)) +
  guides(colour=guide_legend(title=""))
## Zoom: 11

After the Fire

geo_marshall_2=filter(geo_marshall,!is.na(air_quality_2)) %>% 
  st_as_sf(coords=c("lon","lat"), crs = 4326)
marshall.sp.aq2 = as(geo_marshall_2,"Spatial")
factpal <- colorFactor(c("#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641"), c("Strongly disagree","Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"))

leaflet(marshall.sp.aq2) %>%
  addTiles() %>%
  addCircleMarkers(group=marshall.sp.aq2$air_quality_2, color = ~factpal(air_quality_2),
                   radius = 1,
                   opacity = 1,
                   lng = marshall.sp.aq2@coords[,1],
                   lat = marshall.sp.aq2@coords[,2])
ggplot(filter(geo_marshall_2,!is.na(air_quality_2))) +
  annotation_map_tile() +
  annotation_scale() +
  geom_sf(aes(color=as.factor(air_quality_2))) +
  scale_color_manual(values = c("Strongly disagree"="#d7191c",
                                "Somewhat disagree"="#fdae61",
                                "Neither agree nor disagree"="#ffffbf",
                                "Somewhat agree"= "#a6d96a",
                                "Strongly agree" = "#1a9641")) +
  labs(title="Marshall Fire Survey Results (NA's filtered out)", 
       subtitle=paste0("Wave One: Total Respondants Displayed: ",nrow(filter(geo_marshall_2,!is.na(air_quality_2)))),
       caption="Statement: Currently, I am confident that the air in \n my neighborhood is safe to breath") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5)) +
  guides(colour=guide_legend(title=""))
## Zoom: 11

Moran’s I of perception of air quality in one’s neighborhood

library(moranfast)
moranfast(marshall.sp.aq1$air_quality_1, marshall.sp.aq1@coords[,1],marshall.sp.aq1@coords[,2])
## $observed
## [1] 0.01173219
## 
## $expected
## [1] -0.00173913
## 
## $sd
## [1] 0.009821277
## 
## $p.value
## [1] 0.1701736
moranfast(marshall.sp.aq2$air_quality_2, marshall.sp.aq2@coords[,1],marshall.sp.aq2@coords[,2])
## $observed
## [1] 0.003423803
## 
## $expected
## [1] -0.001751313
## 
## $sd
## [1] 0.009905035
## 
## $p.value
## [1] 0.6013408

Map of perception of air quality in one’s home - color from stacked plot colors

Before the Fire

geo_marshall_3=filter(geo_marshall,!is.na(air_quality_3)) %>% 
  st_as_sf(coords=c("lon","lat"), crs = 4326)
marshall.sp.aq3 = as(geo_marshall_3,"Spatial")
factpal <- colorFactor(c("#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641"), c("Strongly disagree","Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"))

leaflet(marshall.sp.aq3) %>%
  addTiles() %>%
  addCircleMarkers(group=marshall.sp.aq3$air_quality_3, color = ~factpal(air_quality_3),
                   radius = 1,
                   opacity = 1,
                   lng = marshall.sp.aq3@coords[,1],
                   lat = marshall.sp.aq3@coords[,2])
ggplot(filter(geo_marshall_3,!is.na(air_quality_3))) +
  annotation_map_tile() +
  annotation_scale() +
  geom_sf(aes(color=as.factor(air_quality_3))) +
  scale_color_manual(values = c("Strongly disagree"="#d7191c",
                                "Somewhat disagree"="#fdae61",
                                "Neither agree nor disagree"="#ffffbf",
                                "Somewhat agree"= "#a6d96a",
                                "Strongly agree" = "#1a9641")) +
  labs(title="Marshall Fire Survey Results (NA's filtered out)", 
       subtitle=paste0("Wave One: Total Respondants Displayed: ",nrow(filter(geo_marshall_3,!is.na(air_quality_3)))),
       caption="Statement: Currently, I am confident that the air in \n my neighborhood is safe to breath") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5)) +
  guides(colour=guide_legend(title=""))
## Zoom: 11

After the Fire

geo_marshall_4=filter(geo_marshall,!is.na(air_quality_4)) %>% 
  st_as_sf(coords=c("lon","lat"), crs = 4326)
marshall.sp.aq4 = as(geo_marshall_4,"Spatial")
factpal <- colorFactor(c("#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641"), c("Strongly disagree","Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"))

leaflet(marshall.sp.aq4) %>%
  addTiles() %>%
  addCircleMarkers(group=marshall.sp.aq3$air_quality_4, color = ~factpal(air_quality_4),
                   radius = 1,
                   opacity = 1,
                   lng = marshall.sp.aq4@coords[,1],
                   lat = marshall.sp.aq4@coords[,2])
ggplot(filter(geo_marshall_4,!is.na(air_quality_4))) +
  annotation_map_tile() +
  annotation_scale() +
  geom_sf(aes(color=as.factor(air_quality_4))) +
  scale_color_manual(values = c("Strongly disagree"="#d7191c",
                                "Somewhat disagree"="#fdae61",
                                "Neither agree nor disagree"="#ffffbf",
                                "Somewhat agree"= "#a6d96a",
                                "Strongly agree" = "#1a9641")) +
  labs(title="Marshall Fire Survey Results (NA's filtered out)", 
       subtitle=paste0("Wave One: Total Respondants Displayed: ",nrow(filter(geo_marshall_4,!is.na(air_quality_4)))),
       caption="Statement: Currently, I am confident that the air in \n my neighborhood is safe to breath") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5)) +
  guides(colour=guide_legend(title=""))
## Zoom: 11

Moran’s I of perception of air quality in one’s home

moranfast(marshall.sp.aq3$air_quality_3, marshall.sp.aq3@coords[,1],marshall.sp.aq3@coords[,2])
## $observed
## [1] 0.01325886
## 
## $expected
## [1] -0.001805054
## 
## $sd
## [1] 0.01040237
## 
## $p.value
## [1] 0.1475826
moranfast(marshall.sp.aq4$air_quality_4, marshall.sp.aq4@coords[,1],marshall.sp.aq4@coords[,2])
## $observed
## [1] 0.009891852
## 
## $expected
## [1] -0.001798561
## 
## $sd
## [1] 0.01042359
## 
## $p.value
## [1] 0.2620606

Map of Wave One symptoms – do a separate map for each type of symptom with a color for yes and a different color for no on that symptom

Map of Wave Two symptoms – do a separate map for each type of symptom with a color for yes and a different color for no on that symptom

After the Fire

geo_marshall_4=filter(geo_marshall,!is.na(air_quality_4)) %>% 
  st_as_sf(coords=c("lon","lat"), crs = 4326)
marshall.sp.aq4 = as(geo_marshall_4,"Spatial")
factpal <- colorFactor(c("#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641"), c("Strongly disagree","Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"))

leaflet(marshall.sp.aq4) %>%
  addTiles() %>%
  addCircleMarkers(group=marshall.sp.aq4$air_quality_4, color = ~factpal(air_quality_4),
                   radius = 1,
                   opacity = 1,
                   lng = marshall.sp.aq4@coords[,1],
                   lat = marshall.sp.aq4@coords[,2])
ggplot(filter(geo_marshall_4,!is.na(air_quality_3))) +
  annotation_map_tile() +
  annotation_scale() +
  geom_sf(aes(color=as.factor(air_quality_3))) +
  scale_color_manual(values = c("Strongly disagree"="#d7191c",
                                "Somewhat disagree"="#fdae61",
                                "Neither agree nor disagree"="#ffffbf",
                                "Somewhat agree"= "#a6d96a",
                                "Strongly agree" = "#1a9641")) +
  labs(title="Marshall Fire Survey Results (NA's filtered out)", 
       subtitle=paste0("Wave One: Total Respondants Displayed: ",nrow(filter(geo_marshall_4,!is.na(air_quality_4)))),
       caption="Statement: Currently, I am confident that the air in \n my neighborhood is safe to breath") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5)) +
  guides(colour=guide_legend(title=""))
## Zoom: 11

Physical Health Symptoms

Wave One:

Wave One counts of symptoms bar plot or pie chart

wave1.symptoms=grep("^symptoms_\\d", names(wave.one), value=T)
rnames=c("Dry Cough","Wet Cough","Wheeze","Itchy or watery eyes",
         "Sore Throat","Headache","Shortness of Breath",
         "Difficult or labored breathing","Sneezing or stuffy nose",
         "Nausea or vomiting ","Allergic skin reaction",
         "Strange taste in your mouth","None of these")
count=list()
for(i in wave1.symptoms) {
  v=sum(!is.na(wave.one[,i]))
  count=c(count,v)
}
df_sym=cbind(wave1.symptoms,count,rnames) %>% 
  as.data.frame() %>% 
  dplyr::select(rnames,count)
df_sym$count=as.integer(df_sym$count)
colnames(df_sym)=c("Symptom","Count")
knitr::kable(df_sym)
Symptom Count
Dry Cough 185
Wet Cough 14
Wheeze 43
Itchy or watery eyes 225
Sore Throat 170
Headache 203
Shortness of Breath 58
Difficult or labored breathing 45
Sneezing or stuffy nose 171
Nausea or vomiting 22
Allergic skin reaction 33
Strange taste in your mouth 75
None of these 387 
# final_df = df %>% 
#   mutate(Percent=round(Count/nrow(wave.one)*100,2))
b=barplot(df_sym$Count, names.arg=df_sym$Symptom, las=2,
          main="Count of Symptoms")
text(b, df_sym$Count-10, df_sym$Count, font=2)

Wave One table of symptom counts by distance category

count.dist=v=wave.one %>% 
  filter(!is.na(.[,wave1.symptoms[1]])) %>% 
  group_by(group) %>%
  count(.[,wave1.symptoms[1]]) %>% 
  as.data.frame() %>% 
  mutate(symptom=rnames[1]) %>% 
  dplyr::select(group,symptom,n)
for(i in 2:length(wave1.symptoms)) {
  v=wave.one %>% 
    filter(!is.na(.[,wave1.symptoms[i]])) %>% 
    group_by(group) %>%
    count(.[,wave1.symptoms[i]]) %>% 
    as.data.frame() %>% 
    mutate(symptom=rnames[i]) %>% 
    dplyr::select(group,symptom,n)
  count.dist=rbind(count.dist,v)
}
colnames(count.dist)=c("Distance","Symptom", "Count")
ggplot(data=count.dist, aes(x=Symptom,y=Count,fill=Distance)) +
  geom_bar(stat="identity") + 
  labs(title="Symptoms")

ggplot(data=count.dist, aes(x=Distance,y=Count,fill=Symptom)) +
  geom_bar(stat="identity") + 
  labs(title="Symptoms")

Wave One table of symptom counts by impact category

count.impact= filter(wave.one,impact_cat!="") %>% 
  filter(!is.na(.[,wave1.symptoms[1]])) %>% 
  group_by(impact_cat) %>%
  count(.[,wave1.symptoms[1]]) %>% 
  as.data.frame() %>% 
  mutate(symptom=rnames[1]) %>% 
  dplyr::select(impact_cat,symptom,n)
for(i in 2:length(wave1.symptoms)) {
  v=filter(wave.one,impact_cat!="") %>% 
    filter(!is.na(.[,wave1.symptoms[i]])) %>% 
    group_by(impact_cat) %>%
    count(.[,wave1.symptoms[i]]) %>% 
    as.data.frame() %>% 
    mutate(symptom=rnames[i]) %>% 
    dplyr::select(impact_cat,symptom,n)
  count.impact=rbind(count.impact,v)
}
colnames(count.impact)=c("Impact_Category","Symptom", "Count")
ggplot(data=count.impact, aes(x=Symptom,y=Count,fill=Impact_Category)) +
  geom_bar(stat="identity") + 
  labs(title="Symptoms - Impact Category",
       subtitle = "I can change the colors of these graphs")

Wave Two:

Wave Two counts of symptoms bar plot or pie chart

wave2.symptoms=grep("^symptoms_\\d", names(wave.two), value=T)[-13:-14] # has 2 more columns
count2=list()
for(i in wave2.symptoms) {
  v=sum(wave.two[,i]==1, na.rm=T)
  count2=c(count2,v)
}

df_sym2=cbind(wave2.symptoms,count2,rnames) %>% 
  as.data.frame() %>% 
  dplyr::select(rnames,count2)
df_sym2$count2=as.integer(df_sym2$count2)
colnames(df_sym2)=c("Symptom","Count")
knitr::kable(df_sym2)
Symptom Count
Dry Cough 88
Wet Cough 12
Wheeze 14
Itchy or watery eyes 93
Sore Throat 53
Headache 64
Shortness of Breath 22
Difficult or labored breathing 4
Sneezing or stuffy nose 81
Nausea or vomiting 8
Allergic skin reaction 18
Strange taste in your mouth 14
None of these 389 
# final_df = df %>% 
#   mutate(Percent=round(Count/nrow(wave.one)*100,2))
b=barplot(df_sym2$Count, names.arg=df_sym2$Symptom, las=2,
          main="Count of Symptoms")
text(b, df_sym2$Count+15, df_sym2$Count, font=2)

Wave Two table of symptom counts by distance category

count2.dist=v=wave.two %>% 
  filter(!is.na(.[,wave2.symptoms[1]])) %>% 
  group_by(group) %>%
  count(.[,wave2.symptoms[1]]) %>% 
  as.data.frame() %>% 
  mutate(symptom=rnames[1]) %>% 
  dplyr::select(group,symptom,n)
for(i in 2:length(wave1.symptoms)) {
  v=wave.two %>% 
    filter(!is.na(.[,wave2.symptoms[i]])) %>% 
    group_by(group) %>%
    count(.[,wave2.symptoms[i]]) %>% 
    as.data.frame() %>% 
    mutate(symptom=rnames[i]) %>% 
    dplyr::select(group,symptom,n)
  count2.dist=rbind(count2.dist,v)
}
colnames(count2.dist)=c("Distance","Symptom", "Count")
ggplot(data=count.dist, aes(x=Symptom,y=Count,fill=Distance)) +
  geom_bar(stat="identity") + 
  labs(title="Wave 2 Symptoms")

ggplot(data=count2.dist, aes(x=Distance,y=Count,fill=Symptom)) +
  geom_bar(stat="identity") + 
  labs(title="Wave 2 Symptoms")

Wave Two table of symptom counts by impact category

count2.impact= filter(wave.two,impact_cat!="") %>% 
  filter(!is.na(.[,wave2.symptoms[1]])) %>% 
  filter(.[,wave2.symptoms[1]]==1) %>% 
  group_by(impact_cat) %>%
  count(.[,wave2.symptoms[1]]) %>% 
  as.data.frame() %>% 
  mutate(symptom=rnames[1]) %>% 
  dplyr::select(impact_cat,symptom,n)

for(i in 2:length(wave2.symptoms)) {
  v=filter(wave.two,impact_cat!="") %>% 
      filter(!is.na(.[,wave2.symptoms[i]])) %>% 
      filter(.[,wave2.symptoms[i]]==1) %>% 
      group_by(impact_cat) %>%
      count(.[,wave2.symptoms[i]]) %>% 
      as.data.frame() %>% 
      mutate(symptom=rnames[i]) %>% 
      dplyr::select(impact_cat,symptom,n)
  count2.impact=rbind(count2.impact,v)
}
colnames(count2.impact)=c("Impact_Category","Symptom", "Count")
ggplot(data=count2.impact, aes(x=Symptom,y=Count,fill=Impact_Category)) +
  geom_bar(stat="identity") + 
  labs(title="Symptoms - Impact Category",
       subtitle = "I can change the colors of these graphs")

Wave One to Wave Two comparison – physical symptoms (maybe select to only be the people who responded to both waves)

compare = left_join(df_sym,df_sym2, by="Symptom")
colnames(compare)=c("Symptom","Wave_One","Wave_Two")
df = melt(compare, id=c("Symptom")) %>% 
  dplyr::mutate(Symptom=rep(1:13,2), 
         Symptom=recode_factor(Symptom,
                               "1"=rnames[1],
                               "2"=rnames[2],
                               "3"=rnames[3],
                               "4"=rnames[4],
                               "5"=rnames[5],
                               "6"=rnames[6],
                               "7"=rnames[7],
                               "8"=rnames[8],
                               "9"=rnames[9],
                               "10"=rnames[10],
                               "11"=rnames[11],
                               "12"=rnames[12],
                               "13"=rnames[13]))
ggplot(df) +
  geom_bar(aes(x = Symptom, y = value, fill = variable), 
           stat="identity", position = "dodge", width = 0.7) +
  scale_fill_manual("",
                    values = c("red","blue"), 
                    labels = c(" Wave One", " Wave Two")) +
  labs(title="Wave One / Wave Two Comparison",
       subtitle="Add counts to top of bars",
       x="Symptom",
       y="Count") +
  theme(plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5),
        plot.caption = element_text(hjust=0.5),
        axis.text.x = element_text(angle = 45, hjust = 1))

Predictors of Physical Health Symptoms

Discussion: FOR LATER